#include "grate_simple.h"
#include "grate_constants.h"
#include <cstring>

double calc_grate_simple(dfa * DFA, const input_data & INPUT, const grate_config & CONFIG){
  int total = DFA->total;

  char * scc = new char [total + 1];
  memset(scc, 0, sizeof(char) * (total + 1));
  
  int scc_num = tarjan_calc_scc(DFA, INPUT, scc);
  int nontriv_scc_num = 0;

  if (CONFIG.OUTPUT_PROGRESS){
    INPUT.dump("The number of non-singleton strongly connected components: %d\n", scc_num);
  }

  double grate = 0;
  for (int cur_scc=1; cur_scc<=scc_num; cur_scc++){
    int scc_size;
    int scc_edgnum;
    double crate = calc_scc(DFA, INPUT, CONFIG, scc, cur_scc, scc_size, scc_edgnum);
    if (crate > grate){
      grate = crate;
    }
    nontriv_scc_num += (scc_edgnum > scc_size);
    if (CONFIG.OUTPUT_PROGRESS){
      INPUT.dump("Growth rate at scc %d: %.20lf\n", cur_scc, crate);
      INPUT.dump("The number of vertices in scc %d: %d\n", cur_scc, scc_size);
      INPUT.dump("The number of edges in scc %d: %d\n", cur_scc, scc_edgnum);
      INPUT.dump("\n");
    }
  }

  if (CONFIG.OUTPUT_PROGRESS){
    INPUT.dump("The number of non-trivial scc: %d\n", nontriv_scc_num);
  }
  delete[] scc;

  return grate;
}

int tarjan_calc_scc(dfa * DFA, const input_data & INPUT, char * scc){
  int total = DFA->total;

  int * mtime = new int [total + 1];
  memset(mtime, 0, sizeof(int) * (total + 1));

  int * father = new int [total + 1];

  int * stack = new int [total + 1];
  int ssize = 0;

  char * nxt_used = new char [total + 1]; // binary encodeing here: 15-mask of number of edges done [+16] = used
  memset(nxt_used, 0, sizeof(char) * (total + 1));

  father[1] = 0;
  mtime[1] = 1;
  int ctime = 1, cver = 1, cedge, nver, pver;
  int scc_num = 0;
  while (cver > 0){
    cedge = nxt_used[cver] & LONG_MASK;
    if (cedge == 0){ 
      nxt_used[cver] ^= 16; // is in stack
      stack[++ssize] = cver;
    }
    else{ // check the son
      pver = (*DFA)[cver][cedge-1];
      if ((pver > 0) && (father[pver] == cver) && (mtime[pver] < mtime[cver])){
        mtime[cver] = mtime[pver];
        nxt_used[cver] |= 32; // improvement
      }
    }
    if (cedge >= INPUT.ALPH_SIZE){
      // handle vertex away
      if ((nxt_used[cver] & 32) == 0){ // no improvement
        if (stack[ssize] == cver){
          nxt_used[cver] ^= 16; // ouf of stack
          ssize--;
        }
        else{
          scc_num++;
          do{
            pver = stack[ssize--];
            nxt_used[pver] ^= 64; // in scc
            nxt_used[pver] ^= 16; // ouf of stack
            scc[pver] = scc_num;
          }while (pver != cver);
        }
      }

      cver = father[cver];
      continue;
    }
    nxt_used[cver]++;
    nver = (*DFA)[cver][cedge];
    if (!nver) continue;
    if (!mtime[nver]){
      ctime++;
      mtime[nver] = ctime;
      father[nver] = cver;
      cver = nver;
      continue;
    }
    else if ((nxt_used[nver] & 16) && (mtime[nver] < mtime[cver])){
      mtime[cver] = mtime[nver];
      nxt_used[cver] |= 32; // improvement
    }
  }
  delete[] mtime;
  delete[] father;
  delete[] stack;
  delete[] nxt_used;
  return scc_num;
}

double calc_scc (dfa * DFA, const input_data & INPUT, const grate_config & CONFIG, char * scc, int scc_id, int & scc_size, int & scc_edgnum){
  int i, j;
  int total = DFA->total;
  scc_edgnum = 0;
  scc_size = 0;

  int * zcc = new int [total + 1];
  memset(zcc, 0, sizeof(int) * (total + 1));

  for (i=1; i<=total; i++){
    if (scc[i] == scc_id){
      scc_size++;
      zcc[i] = scc_size;
    }
  }
  for (i=1; i<=total; i++){
    if (zcc[i] > 0){
      for (j=0; j<INPUT.ALPH_SIZE; j++){
        int v = (*DFA)[i][j];
        if (v > 0 && zcc[v] > 0){
          scc_edgnum++;
        }
      }
    }
  }

  if (scc_size == scc_edgnum){
    return 1.0;
  }

  int * edge1 = new int [scc_edgnum];
  int * edge2 = new int [scc_edgnum];
  int cnum = 0;
  for (i=1; i<=total; i++){
    if (zcc[i] > 0){
      for (j=0; j<INPUT.ALPH_SIZE; j++){
        int v = (*DFA)[i][j];
        if (v > 0 && zcc[v] > 0){
          edge1[cnum] = zcc[i];
          edge2[cnum++] = zcc[v];
        }
      }
    }
  }

  delete[] zcc;

  double * cnt_old = new double [scc_size + 1];
  if (CONFIG.OUTPUT_PROGRESS){
    INPUT.dump("old counters allocated\n");
  }
  double * cnt_new = new double [scc_size + 1];
  if (CONFIG.OUTPUT_PROGRESS){
    INPUT.dump("new counters allocated\n");
  }
  
  for (i=1; i<=scc_size; i++){
    cnt_new[i] = 0;
    cnt_old[i] = 1;
  }
  
  double minrate, maxrate, crate;
  int iter = 0;
  while (iter < CONFIG.MAXITER){
    iter++;
    for (i=0; i<scc_edgnum; i++){
      cnt_new[edge2[i]] += cnt_old[edge1[i]];
    }
    minrate = 1e100;
    maxrate = 0;
    for (i=1; i<=scc_size; i++){
      if (CONFIG.ADD_LOOPS){
        cnt_new[i] += cnt_old[i];
      }
      crate = cnt_new[i] / cnt_old[i];
      if (crate < minrate) minrate = crate;
      if (crate > maxrate) maxrate = crate;
      cnt_old[i] = cnt_new[i];
      cnt_new[i] = 0;
    }
    if (maxrate - minrate < 2 * INPUT.delta){
      double res = (minrate + maxrate) / 2 - (CONFIG.ADD_LOOPS == true);
      if (minrate > MAX_GRATE) res = 0;
      if (CONFIG.OUTPUT_PROGRESS) INPUT.dump("\n%d iterations used\n\n", iter);
      if (GUI_VERSION) printf("%d\n", iter);

      delete[] edge1;
      delete[] edge2;
      delete[] cnt_old;
      delete[] cnt_new;

      return res;
    }
  }
  if (CONFIG.OUTPUT_PROGRESS) INPUT.dump("\nNot enough iterations or precision error! %.20lf precision only!\n\n", (maxrate-minrate)/2);
  if (GUI_VERSION) printf("-1\n");
  double res = (maxrate + minrate) / 2 - (CONFIG.ADD_LOOPS == true);

  delete[] edge1;
  delete[] edge2;
  delete[] cnt_old;
  delete[] cnt_new;

  return res;
}